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QCD at non-zero baryon density is expected to have a critical point where the zero-density cross- 
over turns into a first order phase transition. To identify this point we scan the density-temperature 
space using a canonical ensemble method. For a given temperature, we plot the chemical potential 
as a function of density looking for an "S-shape" as a signal for a first order transition. We carried 
out simulations using Wilson fermions with m % ss lGeV on 6 3 x 4 lattices. As a benchmark, we 
ran four flavors simulations where we observe a clear signal. In the two flavors case we do not 
see any signal for temperatures as low as 0.83T ( . Preliminary results for the three flavor case are 
also presented. 
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1. Introduction 

In recent years, full QCD simulations have become feasible due to development of new al- 
gorithms and increasing computational power. Lattice simulations using dynamical fermions can 
now be performed at finite temperature and zero baryon density. However, simulations at non-zero 
baryon density remain a challenge for lattice QCD due to the complex nature of the fermionic de- 
terminant where the conventional Monte Carlo methods fail. The standard solution of splitting the 
action into the real and positive part and an extra phase fails due to sign and overlap problems. 
To address the overlap problem, a method based on the canonical partition function has been pro- 
posed [0]. While expensive - every update involves the evaluation of the fermionic determinant 
- finite baryon density simulations based on this method proved feasible ^ and a program was 
outlined to scan the QCD phase diagram to look for the critical point [|L ffl]. 

In this paper, we present results based on simulations on 6 3 x 4 lattices with Wilson fermions. 
We plot the chemical potential as a function of baryon density and we clearly observe the "S- 
shape" structure in the Nf = 4 case, indicating a first order phase transition [|p. We do not see such 
a structure in the Nf = 2 case down to 0.8371. We will also present preliminary results for Nf = 3. 

2. Algorithm 

The simplest way to show how to build the canonical ensemble in Lattice QCD is to start from 
the fugacity expansion, 

Z(V,T,^l) = , £Zc(V,T,n)e^ T , (2.1) 

k 

where k is the net number of quarks (number of quarks minus the number of anti-quarks) and Zq 
is the canonical partition function. Using the fugacity expansion, it is easy to see that we can write 
the canonical partition function as a Fourier transform of the grand canonical partition function, 

Z c (V,T,k) = — j o d<Pe-^Z(V,T,n)\^ T - (2.2) 

As an illustration, we will consider the case of two degenerate flavors. After integrating out 
the fermionic part, we get a simple expression 

Zc(V,T,k) = J @\Je- s ^ u) det k M 2 {U), (2.3) 

where 

l r 27t 

det k M 2 (U) = —J d^e-' k UetM(m,n;U) 2 \^T, (2.4) 

is the projected determinant with the fixed net quark number k. de.t k M 2 (U) is a real number due to 
the charge conjugation symmetry of the canonical partition function. However, it is not necessarily 
positive. In our simulations, we use | Redet£M 2 (£/ ) | and fold the phase factor in the observables. 

Exact determinant calculation of fermion matrix is very demanding even on 6 3 x 4 lattices. An 
alternative is to use a noisy estimator [Q, [7]] but this is quite cumbersome. For simplicity sake, we 
used exact evaluation of the determinant in this study. Another technical problem has to do with 
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the Fourier transform; our original approach was to use an approximation where we replaced the 
continuous definition with a discrete one, i.e.: 

det k M 2 (U) « i X>-^MetM(£/ ,) 2 , ^ = ^. (2.5) 

It was shown that the errors introduced by this approximation are small [Q] at least for small quark 
numbers. There are two problems with this approach: the computation time increases linearly with 
the net quark number and, for large enough densities, the Fourier components become too small to 
be evaluated using double precision floating point numbers. To address these issues, in this study 
we used a different approximation method: the winding number expansion This method is 
both faster and more accurate than the original method. 

As mentioned above, we used the chemical potential to probe the phase transition. In the 
canonical approach the chemical potential is measured as the increase in free energy when we 
introduce one more baryon in the system, i.e.: 



_ F(n B + \)-F{n B ) _ 1 Z c (3^ + 3) _ 1 (y(U)) 
WnB (n B + l)-n B ~ /J'" Z c (3n B ) ~ j3 (a(U)) 



where 



, x Redet 3;lR M 2 ([/) 

a(U) = 8 .) , and (2.7) 

V ' |Redet 3 , lfl M 2 (£/)| 

Redet 3 „ B+3 M 2 (t/) 
Y[U) |Redet 3 „ fl M2(f/)| ' ^ 

a(U) is the phase and () stands for the average over the ensemble generated with measure 
|Redet 3 „ 6 M 2 (£/)|. If the configurations in the ensemble have equal probability for both positive 
and negative a(U), the denominator in the equation above becomes zero and we have a "sign 
problem". In our simulations the sign oscillations are under control. 



3. QCD phase diagram 

At zero baryon density, it has been known for quite some time that QCD undergoes a transition 
from a confined phase to a deconfined phase at a temperature T c « 170MeV. Lattice QCD suggests 
that the transition is in fact a smooth crossover. This is expected to turn into a first order phase 
transition as the baryon density is increased. A schematic picture of the expected phase diagram in 
the density-temperature plane (see Fig. [|) shows the crossover ending with a critical point at some 
non-zero baryon density. The first order phase transition is characterized by a coexistence region 
separating the hadronic phase and the plasma phase. 

To search for the phase boundaries of the coexistence region, we scan the phase diagram by 
varying the density while keeping the temperature fixed. The baryon chemical potential should 
exhibit an "S-shape" as one crosses the coexistence region [Eft. 

In our simulations we use quarks that are much heavier than the physical quarks. Furthermore, 
for simplicity sake, we carry out simulations where the quark masses are degenerate. We will 
present results for simulations using 2, 3 and 4 degenerate flavors of quarks. The expected phase 
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Figure 1: Schematic phase diagram of QCD. 




Figure 2: Expected phase diagrams of two flavors (left) and four flavors (right). 



diagrams for 2 and 4 flavor cases are shown in Fig. ^[ The two flavor case is expected to have a 
diagram very similar to the full QCD one whereas it is known that in the case of four flavors the 
first order phase transition extends all the way to zero baryon density. We will use the four flavor 
simulations as a benchmark to show that the methodology we use can determine the boundaries of 
the coexistence region. 

4. Results 

In Fig. H we show the results for our four flavor simulation. On the technical side, we note 
that we didn't have a sign problem: even the simulations at the largest density where the box size 
is 1.8fm and the temperature is 0.90r o the sign oscillations were moderate. The plots show a clear 
signal for a first order phase transition when the temperatures are lower than T c . 

To identify the boundaries of the coexistence region and the critical value for the chemical 
potential we used the Maxwell construction More precisely, we selected four points in the 
"S-shape" region and we fit these points with a third order polynomial. A better approach would 
be to use some phenomenologically motivated functional form and try to fit a larger region; we 
changed the fit function and we also extended the fit region. For all reasonable fits, we found that 
the values of the boundary points, pi and P2, and the value of the critical chemical potential, pL c , are 
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Nf=4, Wilson gauge and fermion action, 6 x4 




2 4 6 8 10 

Baryon number 



Figure 3: Baryon chemical potential vs. baryon number at different temperatures for Nf = 4. 



fairly insensitive to our choice of the fit function or fit region - the simple third order polynomial 
fit was sufficient. 
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Figure 4: Maxwell construction for T = 0.942" e and Nf = 4. 



Once the pi and P2 are determined for different temperatures, they can be used to plot the 
boundaries of the coexistence region. We can also find the critical point by determining where the 
width of the coexistence region shrinks to zero. In Fig. |5| we compare our results to those from 
a study using staggered fermions [g]. We see that our coexistence region is much narrower. This 
could be due to our heavier pion mass (m K m lGeV compare to m K 300MeV), or due to the fact 
that we use a different fermion formulation. 

The fact that the results of our four flavors simulations are consistent with the expectation 
and with other lattice studies is encouraging. The only issue that needs to be addressed is the 
discrepancy in the location of the boundaries. 

In Fig. |6] we show our results for Nf = 2 simulations. These simulations are more expensive 
than the four flavors simulations due to sign fluctuations. In the two simulations that we ran at 
T = 0.86!T C and at T = 0.83T f we do not see any signal for a first order phase transition. There is at 
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Figure 5: Phase boundary for four flavors QCD: our work and the results of a study using staggered 
fermions [B]. (Dashed lines are plotted to guide the eye; they are not the result of an extrapolation. The 
error bars are not determined yet.) 

least one claim [Q] that the critical point occurs at temperatures as low as T = O.ST c . If this is indeed 
the case, we need to run simulations at even lower temperatures in order to see the "S-shape". 



Nf=2, Wilson gauge and fermion action, 6 3 x4 
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Figure 6: Baryon chemical potential vs. baryon number for Nf = 2. 

Three quark flavors are relevant for the structure of matter at energies of the order of few 
hundred MeV: two light quarks u, d and one heavier quark s. Lattice simulations that come close 
to approximate full QCD treat the lighter quark flavors as degenerate and use one heavier quark 
flavor. Ideally, we would like to carry out simulations close to the physical point. However, this is 
not really practical especially for Wilson fermions on lattices as coarse as we use in our study. In 
fact, our quark masses are even heavier than the strange quark mass. We decided to investigate the 
Nf = 3 with the hope that the phase diagram is, at least qualitatively, close to the full QCD phase 
diagram. 

In Fig. [7| we show the results of our simulations for Nf = 3. For these simulations we used 
Iwasaki gauge action and clover fermions in order to reduce lattice discretization errors. For T = 
0.92r c , we don't see any signal for a first order phase transition. 
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Nt=3, Iwasaki + clover action, 6 3 x4 
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Figure 7: Baryon chemical potential vs. baryon number at different temperatures for four flavors. Simula- 
tions are done using Iwasaki gauge action and clover fermion action. 

5. Summary 

We presented results for QCD simulations at non-zero baryon density using the canonical 
ensemble approach. In the four flavor case we see a clear signal for a first order phase transition - 
this is consistent with phenomenological expectations and previous lattice results. This proves that 
our method is sound. In the two flavor case, we do not see any signal for temperatures as low as 
0.83r c . We also presented results for three flavors where we do not find any signal for a transition 
at T = Q.92T C . We plan to continue our investigations at lower temperatures and smaller quark 
masses. 
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